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ABSTRACT 



We employ the PLUTO code for computational astrophysics to assess and compare the validity of different numerical algorithms on 
simulations of the magneto-rotational instability in 3D accretion disks. In particular we stress on the importance of using a consistent 
upwind reconstruction of the electro-motive force (EMF) when using the constrained transport (CT) method to avoid the onset of 
numerical instabilities. We show that the electro-motive force (EMF) reconstruction in the classical constrained transport (CT) method 
for Godunov schemes drives a numerical instability. The well-studied linear growth of magneto-rotational instability (MRI) is used 
as a benchmark for an inter-code comparison of PLUTO and ZeusMP. We reproduce the analytical results for linear MRI growth in 
3D global MHD simulations and present a robust and accurate Godunov code which can be used for 3D accretion disk simulations in 
curvilinear coordinate systems. 
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1. Introduction 

Most classical work on the MRI has bee n conducted us- 
ing local shearing sheet simulations (Brandenburg et all 119951 : 
iHawlev et all 119951: ISano et dll2004 . Nevertheless, recent de- 
velopment has shown that critical questions about turbulence 
in accretion disks can only be addressed in global simula- 
tions, e.g. the saturation level of turbulence, the global evolu- 
tion and the role of the global field configuration. So far most 
existing codes which can handle global disk simulations are 
based on a Zeus like finite difference scheme (IHawlevl l2000t 
iFromang & Nelson , 2006, 2009). Thus, it would be very bene- 
ficial to have an independent method to attack the open ques- 
tions in Magneto-Hydrodynamic disk sim ulations. The recently 
introduced MHD Godunov code PLUTO (iMignone et al.L 120071) 
brought new insight into the field of jet propagation and for 
instance MHD simulations in the solar system, e.g. the space 
weather prediction. The latter example has the benefit that the 
predictions obtained with the numerical method can be tested 
against observation. A number of investigators have recog- 
nized the importance of using conservative Godunov-type up- 
wind schemes rather than non-conservative finite difference al- 
gorithms which do not make use of the characteristic of the as- 
sociated Riemann Problem. Clearly accretion disk calculations 
using conservative schemes are of great interest. As we want 
to apply our simulations to study the physical conditions for 
the planet formation process in circumstellar disks, we need 
not only the dynamical evolution of the disk, but also the de- 
tailed thermal and chemical conditions of the disk gas. None 
of the global disk simulations so far have been performed with 
tracing the temperature changes in the gaseous disk. The first 
result in this direction is comparing MRI linear growth rates 
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including ra diation tran s port w hich was achieved in local-box 
simulations ( iFlaig et all |2009|) . It is widely accepted that the 
magneto-rotational instability (MRI) drives the turbulence in 
proto-planetary disks, leading to gas accretion onto the star as 
well as radial and vertical mixing. The efficiency of MRI-driven 
turbulence depends on the ionization degree of the gas. To de- 
termine the latter, thermal evolution including heat input from 
magnetic field dissipation has to be taken into account. Therefore 
we need an accurate conservative numerical scheme for treating 
such MHD processes as magnetic diffusion and reconnection in 
the presence of the small-scale turbulence. One of the most im- 
portant challenges is an accurate and effective way to control 
and monitor the constraint of the divergence-free magnetic field. 
In the last years, several approaches on the numerical treatment 
of div(B) have been proposed. Two main classes for Godunov 
schemes have been established over the years. First, Powell's 
"Eight Waves" method uses the modified Euler equations with 
specific source terms and allows for magnetic monopoles treated 
as the 8* wave (Pow ell 1994). This idea i s similar to the ellipti- 
cal cleaning method ( iDedner et al.L |2002|) which damps numer- 
ical monopoles. The second class is the Constrained Transport 
method, initially u sed as MocCT in ZEUS and then adapted to 
Godunov schemes dBalsara & Spiceiill999l) . Here, the magnetic 
field components are always reprojected on the staggered grid 
and div(B) = is kept within machi ne accuracy. The classi- 
cal CT method for Godunov schemes dBalsara & Spiceii 119991) 
uses an arithmetic averaging of the electro-motive forces (EMF). 
This average proves to be not consist ent with the plane par- 
allel pr opagation in upwind scheme s dLondrillo & Del Zannal 
I2OOOI) . dLondrillo & del Zannal |2004|) have proposed a new way, 
named upwind constrained transport (UCT), of combining CT 
and a consistent upwind EMF reconstruction. [Gardiner & Stond 
d2005l) also developed a consistent method for the CT implemen- 
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tation in the ATHENA code. The PLUTO code o ffers a choice of 
the m ethods described above. The PLUTO code (iMignone et al.L 
l2007b is a conservative multi-dimensional and multi-geometry 
code, with applications for Newtonian, relativistic, MHD and 
relativistic MHD fluids. The possibility of switching between 
several numerical modules in PLUTO allows us to perform a 
broad variety of test problems and choose the module with the 
best convergence and stability. E.g. the code offers the possibility 
of using Riemann solvers of ROE, HELD, HEEC, HEE and Lax- 
Friedrich Rusanov. In this paper we investigate the applicability 
of the MHD methods for 3D curvilinear global simulations of the 
turbulent disk. To have comparable and reliable results, simula- 
tions with the ZEUS code are included in this paper, representing 
the group of the finite difference codes. ZeusMP version lb has 
been used successfully fo r 3D MHD global s imula t ions of astro- 
physi cal disks. ZeusMP dStone & NormanL Il992t iHayes et all 
I2OO6I) has a classical staggered discretization of magnetic and 
electric vector fields (MoC CT). An identical model is set up 
in the actual ZeusMP code as well as the PLUTO 3. We ex- 
amine the performance of the different numerical schemes by 
comparing the growth of the MRI with the expectations from 
linear analysis. The PLUTO code has successfully reproduced a 
series of standard MHD tests, including t he rotating shock-tube , 
the fast rotor and the blast wave solution (iMignone et al .1120071) . 
In addition, a magnetized accretion torus in 2.5 dimensions is 
demonstrated. 

We are testing the PLUTO code specifically for the linear MRI. 
Our 3D disk setup allows to reproduce exactly the conditions 
for MRI typical for accretion disk studies. At the same time, 
unphysical results can easily be discovered by analytical linear 
MRI analysis. In section 2 we describe the differences in the nu- 
merical schemes of ZEUS and PLUTO. Section 3 presents our 
disk model used for linear MRI simulation and the measurement 
method. In Section 4 we show the results of the inter-code com- 
parison. A convergence check for the chosen code configuration 
is investigated in section 5. Conclusion and outlook are given in 
section 6. 



2. Numerics 

2.1. MHD Equations 

The equations of ideal magnetohydrodynamics written for 
Godunov schemes in conservative form are: 



— + V ■ (puu - BB) 

dt ^ 



-VP* + pW, 



dp 
dt 



-H V ■ (pH) = 0, 



^ + V iuB-Bu) = 0, 
dt 



and 

dE 

— + V ■ ((£ + P*)u 
dt 



B{Bu)) = puV*!' 



(1) 



(2) 



(3) 



(4) 



where p is the gas density, pu is the momentum density, B is the 
magnetic field, and E is the total energy density. The total pres- 
sure is a sum of magnetic and gas pressure, P* - P + (B ■ B)/2. 
Total energy density is connected to internal energy e as £ = 
e + pu- 12 + 12. Gravitational potential ^ = 1//? is chosen in 
the normalized form to describe the sheared azimuthal flow of 



the gas, reproducing the Kepler rotation. We use cylindric co- 
ordinates in our models with the notation: {R,(p,Z). The equa- 
tions are solved using uniform grids. In our runs we take an adi- 
abatic approximation for the Equation of State, P - {y - l)e, 
with y = 5/3. 

Treatment in ZEUS: The equations (1-4) are solved by ZeusMP 
in a non-conservative way, i.e. the momentum conservation will 
be treated in the following way. 



p— = -VP + pV^ + —[V xB]xB. 
dt An 



(5) 



As the vector components of velocity and magnetic fields are 
stored at the grid interfaces, fluxes can be calculate directly in 
ZEUS. Density and pressure are stored at the cell center and the 
vector values are at the interfaces. The scheme is second order 
in space and time. For time integration, ZEUS uses a Hancock 
scheme. Contrary to Zeus-type finite difference schemes, a 
Godunov type scheme, such as represented by the PLUTO code, 
follows a conservative model. Therefore the temperature can be 
exactly calculated from the total energy conservation. In gen- 
eral all quantities are stored at the cell center and an 2nd order 
Runge-Kutta (RK) iterator is employed. Prior to flux calcula- 
tion, variables are reconstructed at the grid interface, which im- 
plies interpolation with a limited slope. The resulting two differ- 
ent states at the interfaces allow to solve the Riemann problem, 
which is the main feature of a Godunov code. CFL value is 0.25. 

2.2. Magnetic fields in Godunov scliemes 

In case of MHD, cell-centered Godunov schemes have in gen- 
eral problems with occu rring numerical magnetic monopoles 
( lBalsara&Spiceilll999l) . Powell's "8 wave" method includes 
this numerical magnetic monopoles as source terms. In the con- 
strained transport (CT) MHD method, the magnetic fields are 
located at the interfaces. Usually, this method allows to handle 
the magnetic monopoles at machine accuracy. After updating the 
staggered magnetic field, the cell centered variables will be inter- 
polated from the new interface values. The mainly used config- 
uration in curvilinear coordinates systems consists of the second 
order upwind scheme with piece-wise linear reconstruction and 
RK 2 for time evolution. In our inter-code comparison we test 
Powell's "8 Waves" method, t he CT scheme and the different 
MHD Riemann solv ers of R oe (ICargo & GallicellI997h . HELD 
(Mivoshi & KusancJ l2005l) , HEE and Lax-Friedrich Rusanov. 
We must note that for the most frequently used configurations 
the ZEUS code has a speed up of factor 4. This factor appears 
because of the additional time step in the Runge Kutta time inte- 
gration and the solving of the Riemann problem. 

2.3. CT EMF reconstruction 

The solenoidality of magnetic fields is difficult to represent cor- 
rectly in Godunov schemes. The simple arithmetic average to re- 
construct the EMF is presen t ed in F ig. 1 and follows the Eq.(7)- 
Eq.(9) in iBalsara & SpiceJ (Il999l). We share the opinion pre- 
sented in lEondrillo & Del Zannal(f2000h : lLondrillo & del Zannal 
I2OOI), that such a simple hybrid method as ACT in Godunov 
scheme is not solving the monopole problem. The arithmetical 
CT algorithm is not consistent with the integration algorithm 
for plane-parallel, grid-aligned flows, as in our case with the 
dominating Kepler rotation. In PL UTO, t he upwind CT scheme 
(UCT) follows Gardiner & Ston^ (l2005l) and is based on the 
work of iLondrillo & Del Zanna "(120001) : ILondrillo & del Zannal 
(I2OOI . The main difference with respect to ACT is to construct 
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Fig. 1. The reconstruction of the EMF (green axis) from four 
Godunov fluxes (red arrows) in case of arithmetic average of 
fluxes (ACT) 



the EMF update i n such away that it is reco vering the proper di- 
rectional biasing (iGardiner & Stonel i2005i) . The first method is 
here called UCT. 

^z,/+l/2,j+l/2 = 2'-^Z''+l/2j ^z,i+l/2,j+\ + Sz,i,j+l/2 + Sz,MJ+l/2) 
1 

The term i/2,j+ 1/2 is taken from the Godunov fluxes as illus- 
trated in Fig. 1. e~jj are the finite volume electric fields calcu- 
lated from the cell center values. In addition to this method, we 
also test the EMF corner-integration, here called UCTcontact: 

fz,/+l/2,;+l/2 = -(ec,i+l/2,j + &,/+l/2,j+l + e-.ij+l/l + ez.i+l.j+i/l) 
+ ^((-^)!+l/2j+l/4 - (-^)/+I/2,;+3/4) 



6x , , 



with 



+ Y(("^)/+l/4.i+l/2 - (■^)i+3/4J+l/2) 

de^ 2 



and 



\dylM/2.j+l/4 



' (deJdy)ij+i/4 
+(deJdy)i+ij+i/4), otherwise 



iorVjcj+i/2,j >0,\ 
for v^-,+i/2j < 0, 



The latter shows in general the same MRI evolution and therefor 
we concentrate on testing of UCT. 

3. Linear MRI in giobai disl< 



The magneto-rotational instability (MRI) dBalbus & Hawlevl 
[T99Tt iHawlev & BalbusL [Mil iBalbus & Hawlevl [19981) is be- 
lieved to occur in astrophysical disks. The analytical descrip- 
tion of its linear stage has been discovered bv lBalbus & Hawlevl 
(II99 lb . Local shearing box simulations have confirmed the ana- 
lytica l expectation o f the lin ear MRI growth jBrandenburg et all , 
1 19951: iHawley et all Il995h . Global simulati ons and also non- 
linear evolution of the MRI was presented in lHawley & Balb uj 
(Il99lh . Our simplified setup here is representative for a global 
proto-planetary disk and made for a modest radial domain. 



3.1. Global setup 

Our disk model is not representing the whole 3D accretion disk, 
but a part which is large enough to fulfill global properties. The 
gas is rotating differentially with the Keplerian shear, - j. 
There is no vertical stratification in both rotation and gas den- 
sity. Initial temperature and density are constant in the entire 
disk patch, p = po, T - Tq, cq - 0.1. Therefore the or- 
bital speed at the inner boundary in units of the sound speed 
is 10. The domain extends pi/3 radians in the azimuthal direc- 
tion, Z = +0.5 and from 1 to 4 units in the radial direction. 
An annulus of uniform vertical magnetic field is placed at radii 
between 2 and 3 units, avoiding the radial boundaries so that 
any effects of the boundary conditions are reduced. The resolu- 
tion is [R,(f),Z] = [128,64,64]. For both codes we use identi- 
cal random generated velocity perturbations (lO^^WKep) for ra- 
dial and vertical velocities as the initial seed for the turbulence. 
Boundary conditions are periodic for the vertical and azimuthal 
direction, and zero gradient for the radial one. Using the ana- 
lytical prediction for the critical MRI wavelength, we choose 
the uniform vertical magnetic field strength in order to obtain 
the number of fastest-growing wavelengths fitting in the domain 
he igh, = 0.0551 3 / w whe re n = 1,2, ...8 (after equation 2.32 
in iBalbus & Hawlevl (Il99lb ). The critical wavelength and the 
growth rate can be directly measured from our simulations and 
compared with the analytical prediction. 



3.2. Limits of the MRI growth rate 

Before using the growth rate as numerical benchmark, the phys- 
ical processes limiting the growth rate in the linear MRI have 
to be considered. When magneto-rotational instability sets in, 
the magnetic fields are amplified with B = Bq exp(yr) until they 
reach the saturation. Here we use a standard notation where y 
is the MRI growth rate. As it is mentioned in lHawley & BalbusI 
(Il991h . the absolute limit of the growth rate for ideal MHD is 
given in case of zero radial wave vectors qr = with the normal- 
ized wave vector q~. 



q- = A:-. V16/15va/Q. 



(6) 



with the Alven speed va - B/ -\j4-7Tp. Then, the critical mode 
q^ = 0.97 grows with 7 = 0.75Q, visible in Fig. 2. Like in 
Balbus & Hawley (1991b) we do not a priori know the radial 
wavenumber in our simulations. In Fig. 2 we plotted the growth 
rate for the amplitude of different vertical modes. The simulation 
is made for B~ = 0.05513/4 with HELD solver for resolution 
128x64x64 (see also Table 2, HLEDlin)- We have obtained the 
MRI g rowth rates which are comparable with^ Hawlev & BalbusI 
( ll99ll) .(Fig. 8). The addition limi ters of the MRI growth rate are 
viscous and resistive dissipation dPessah & Chanll2008l) . We in- 
terpret the differences between the results, such as number of 
resolved MRI modes and the growth rate, as the effect of numer- 
ical dissipation which is individual for each solver. In case of 
MRI, the numerical dissipation is the reason why we need some 
minimal am ount of the gr i d cells to resolve the fastest growing 
MRI mode. iHawlev et al.l ( Il995h claim to need at least 5 grid 
cells per fastest growing mode. Generally, it holds that the more 
grid cells we use per critical mode the less is the numerical dissi- 
pation. But beside increasing resolution, changing the Riemann 
solver or the reconstruction method also affects the dissipation. 
The piece-wise linear method to reconstruct the interface states 
for the Riemann solver drives in general to higher numerical dis- 
sipation because of the additional interpolation. In ZEUS the ve- 
locity and magnetic fields are akeady at the interface and the 
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Fig. 2. Analytical growth rate is plotted for different q,- (solid 
to dash-doted lines). The triangles present the measured MRI 
growth rates from our simulation with initial axisymmetric white 
noise disturbance from Model HLLDlin 128x64x64 with n=4. 
The values are measu red for each radii and aver aged. The results 
are similar to those in lHawlev & Balbusl(ll991h .(Fig. 8) 
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0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 
Orbits 

Fig. 3. Maximal amplitude growth of the critical MRI mode is 
plotted for the 42 radial slices between radii R - 2 and R - i 
(solid lines). The setup is constructed to excite the mode with 
/lent = H {n - 1). For the radii in the middle of MRI-active do- 
main, the MRI growth rate reaches the value about y - Q.15Q.. 
At the outer edge of the radial domain, the growth rate becomes 
higher than 0.75f2. We suggest that the MRI modes there are 
affected from the inner part where MRI already breaks into non- 
linear evolution. 




2 4 6 8 10 12 

Orbits 



Fig. 4. The total domain integrated radial magnetic energy plot- 
ted for n=4 (top) and n=8 (bottom) with upwind CT (solid line), 
Powell 8 wave (dashed line) and arithmetic CT (dotted line). 
The values are normalized over the initial orbital kinetic energy. 
After the initial disturbance, the growth of critical mode dom- 
inates and there is a linear MRI region visible between 2 and 
8 inner orbits. For HLLD and Roe, the arithmetic CT shows a 
growth stronger then MRI until 5 inner orbits caused by a nu- 
merical instability. HLL and TVDLF have a much lower growth 
because both methods are viscous and do not resolve the critical 
mode for the adopted resolution. 



fluxes can be calculated directly. For the inter-code comparison 
we use the setup described in section 3.1 with a vertical mag- 
netic field with n - 4. Then the resolution in Z is 16 grid cells 
per fastest growing mode which is also enough for very diffuse 
solvers. For better discerning between the code configurations 
we include also a weaker vertical magnetic field with « = 8. 

3.3. Measurement method 

To reach the ideal limit and to explain our measurement method 
we set up a model with 1 critical mode (64 Grid cells per mode) 
to reduce the numerical dissipation (Fig. 3). In addition we use 



here a critical mode initial disturbance to enforce the MRI fastest 
mode to grow faster then any other waves. Fig. 3 shows the am- 
plitude of the critical mode A- visible in the radial magnetic 
field Br{Z). For each radial slice we use the local orbital time 
Tr - IttR^-^ and calculate the growth rate between 1 and 1.5 lo- 
cal orbits. Most radial slices plotted in Fig. 3 show a growth rate 
close to the analytical limit. The inner annuli reach non-linear 
amplitudes first, producing disturbances that propagate outward 
and affect the development of the slower-growing linear modes 
in annuli further from the star. As a result the annuli near R=3 
grow faster than indicated by the linear analysis in all the calcu- 
lations with sufficient spatial resolution. 
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Solver 






8 IV, 


=4 








8W„=8 


ZEUS 


0.69 ± 0.03 








0.69 ± 0.03 








HELD 


0.68 ±0.17 


0.68 ± 0.03 






0.86 ±0.18 


0.70 ± 0.03 






ROE 


0.70 ±0.18 


0.68 ± 0.03 


0.68 ± 


0.02 


0.84 ±0.18 


0.70 ± 0.03 


0.70 


±0.04 


HLL 


0.64 ± 0.05 


0.63 ± 0.05 


0.61 ± 


0.05 


0.55 ± 0.09 


0.51 ±0.06 


0.50 


±0.06 


TVDLF 


0.64 ± 0.05 


0.63 ± 0.05 


0.63 ± 


0.05 


0.55 ± 0.09 


0.51 ±0.06 


0.50 


±0.06 



Table 1. MRI growth rates are measured for models with 4 and 8 MRI modes (n - 4, 8). Measurement of growth rates follows the 
description in section 3.3, y is averaged from the different radial shells. The growth rates higher than the analytical value are marked 
with red. The uncertainty numbers are measured from the radial variations in the growth rate visible in Fig. 7 or Fig. 9. 



4. Inter code comparison 

4. 1. Energy evolution 

Table 1 demonstrates the collection of models we made for inter- 
code comparison. The linear regime of MRI is visible in the to- 
tal radial magnetic energy evolution. Fig. 4. The energy growth 
due to the excitation of the critical mode dominates after 2 in- 
ner orbits and the linear growth regime becomes visible. Fig. 4 
shows the radial magnetic energy evolution for the 4 and 8 mode 
runs. In the 8 mode setups we can clearly distinguish 3 differ- 
ent classes of evolution. The Lax-Friedrich and HLL Riemann 
solver with all possible method of EMF reconstruction (ACT or 
CT with arithmetic EMF averaging, upwind CT and "8 Wave") 
have the slowest energy growth. The second group is built with 
the HLLD and Roe solvers with the upwind CT and the 8 wave 
method. Their results are most similar to ZEUS. The third group 
are the Roe and HLLD solver in combination with ACT (red col- 
ors in Table 1). The steep slope indicates a numerical problem. 
Table [T] presents the growth rate values for the different code 
configuration. The growth rates obtained from the model using 
the ACT scheme are higher than the analytical value. 

4.2. Arithmetic CT 

In case of ACT and the high-resolution solver, the strong growth 
of magnetic energy overshoots the analytical expectation for 
MRI growth (Table 1 and Fig. 4). Here a numerical instabil- 
ity dominate s in the cases with the HLLD and Roe solvers. 
The fastest-growing mode has a vertical wavelength shorter 
than in the linear analysis. In Fig. 5 we plotted the growth 
rate for each radial slice, calculated from I to 1.5 local or- 
bits. Roe and HLLD show the same unphysical high growth 
rate. On first inspection, the more diffusive HLL and TVDLF 
solvers seem not to be affected. However, at higher resolu- 
tions the numerical instability occurs there as well. The in- 
stability shows a 'checkerboard' pattern (Fig. 5). The occur- 
rence of such instabilities was already mentioned in Miyoshi 
& Kusano (2008). The Roe and HLLD solvers which include 
the Alfven characteristic accurately evolve the disturbances in- 
dependently of the resolution. The reason of the numerical insta- 
bil ity Ues in the inconsistent EM F reconstruction, as described 
in iLondrillo & Del Zanna' ('2000'). The correction meth ods are 
propos ed i n [Londrillo & de l Zanna (2004), Gardiner & Stond 
(l2005h and IZieglerl (l2004h . In short one can summarize, that 
the simple arithmetic average of 4 Godunov fluxes to recon- 
struct the EMF is not consistent with the upwind scheme. The 
error becomes visible in case of dominating flows along the az- 
imut hal direction in 3D, pres ented in the Loop advection test 
from iGardiner & Stond (l2005l) . It is interesting to note that the 
results of HLLD and Roe are very close. Also the values of HLL 



ACT, 4 MRI modes 




2.4 2.6 
Radius [AU] 

ACT, 8 MRI modes 




2.4 2.6 
Radius [AU] 

Fig. 5. Local growth rate for n=4 (top) and n=8 (bottom) with 
arithmetic CT. For HLLD and Roe the local growth rate exceeds 
the analytical limit. The problem does not affect the HLL and 
TVDLF solver for the used resolution 



and TVDLF are exact the same in our tests. Resolving the Alfven 
characteristic in MHD Riemann solver plays an important role 
for numerical dissipation. 



4.3. Upwind CT vs 8 wave 

With the Upwind CT EMF reconstruction the critical MRI mode 
dominates and 'checkerboard' pattern disappears. Fig. 6 presents 
a contour plot of the azimuthal magnetic field. HLLD with up- 
wind CT and ZEUS show 4 evolving modes. Fig. 7 shows the 
small difference between "8 Waves" and upwind CT. The ad- 
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Fig. 6. Contour plot of azimuthal magnetic field for arithmetic 
CT (top), UCT (middle) and ZEUS (bottom). In case of arith- 
metic CT, the numerical instability dominates MRI. In case for 
UCT and ZEUS the n=4 is visible. A change from blue to red 
means a change of magnetic field sign. The snap-shot of mag- 
netic field is taken after 8 inner orbits. 



ditional EMF corner-integration UCTcontact does not change 
the result. The strongest effect on numerical dissipation has the 
change of Riemann solver. The reason is the number of wave 
characteristics which is used by each solver. HLL and TVDLF 
use the fast magneto-sonic wave characteristic as maximal signal 
velocity. HLLD includes additional the Alfven wave character- 
istic and the contact discontinuity. The ROE solver includes all 
7 MHD waves, as well as the slow magneto-sonic characteristic. 
ZEUS shows the MRI development same as HLLD and Roe. 

4.4. 4 to 8 MRI modes 

In case of 4 modes, the resolution in the Z direction is suffi- 
cient for all Riemann solvers to resolve the critical mode. HLLD 
, Roe and ZEUS shows the highest growth rate followed by 
TVDLF and HLL. As already presented in Fig. 4. The models 
for 8 MRI critical modes are summarized in Table L Energy 
evolution is shown in Fig. 4 (bottom). Growth rates are signif- 



Fig. 7. Local growth rate for n=4 (top) and n=8 (bottom) with 
upwind CT (solid line), UCTcontact (dashed line) and Powell's 
Eight Wave (dotted line). The high resolution MHD Riemann 
Solver HLLD and Roe show the same MRI evolution, producing 
same growth rate as ZEUS. HLL and TVDLF are more diffusive. 



icantly reduced only for HLL and TVDLF solvers. ZEUS and 
the MHD Riemann solver ROE and HLLD resolve 7 modes, 
HLL and TVDLF only 4 modes. The resolution about 8 grid cells 
per critical mode is not enough to resolve the modes perfectly. 
Nevertheless, the additional Alfven characteristic included in the 
HLLD solver drives to lower numerical diffusion. Fig. 8 shows 
a contour plot of radial magnetic field after 7 inner orbits for 
HLLD solver and ZEUS. The ZEUS code shows the MRI pat- 
tern very close to those of HLLD and Roe. 

4.5. PPM vs PLM 

The reconstruction of the states at the interfaces leads to the high 
dissipation in the Godunov schemes (see section 2). Therefore 
we test the fourth order piece-wise parabolic method (PPM) by 
Collela (1984) for the ROE,HLLD and for the diffusive TVDLF 
solver. Fig. 9 shows the difference in local MRI growth rate 
for simulations made for interpolations with piece-wise linear 
method (PLM) and PPM. The results of MHD Riemann solver 
HLLD and ROE are not strongly affected, though there is a 
slightly lower numerical dissipation visible when PPM is used. 
But there is a strong effect of the interpolation method visible 
for the classical Lax-Friedrich and HLL solvers. In case of the 
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Fig. 8. Contour plot of radial magnetic field in code units for the 
HLLD solver (top) and ZEUS (Bottom). After 7 inner orbits the 
inner radial part breaks into nonlinear regime. This affects also 
the growth rate at the outer part as we have already seen. 



8 mode run, also the classical Riemann solver presents 7 modes 
and the values of ROE, HLLD and TVDLF converge. To verify 
this conclusion we investigate the resolution convergence, in- 
cluding only the HLLD and TVDLF solver with upwind CT and 
both reconstruction methods. 



5. Convergence test 

The models for the convergence tests are summarized in Table 
2. We take here the n=4 model. We double the resolution up 
to 4 times for each dimension to check for the earliest conver- 
gence. In these runs we use an initial axisymmetric white noise 
for the radial and vertical velocity {www.random.org). This leads 
to slightly different growth rate values as in Table I. Fig. 10 
shows the results. The HLLD MHD Riemann solver in combina- 
tion with piece-wise parabolic reconstruction shows the earliest 
convergence. Even with a resolution about 32 in Z and with the 
PPM reconstruction, we can resolve the critical mode and get the 
converged value of about y - 0.700 (Table 2). The resolution 
runs marked with a star do not resolve the fastest growing MRI 
mode. As expected, the Lax Friedrich solver shows the highest 
numerical dissipation. Here 8 grid cells per fastest growing mode 
are not enough. For the highest resolution of 512x256x256 we 
see in all cases a decreasing of the growth rate. We cannot say 
whether this behavior is connected with the presence of high ra- 
dial wavenumber or a change in the numerical dissipation. All 
except the Zeus runs decline by less than the uncertainty. 
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Fig. 9. MRI growth rate for PPM (dashed) and PLM (solid). In 
the case of the Lax-Friedrich solver, the PPM interpolation leads 
to a lower numerical dissipation and growth rates nearer the an- 
alytic solution. Differences between HLLD and ROE solvers are 
negligible. 
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Fig. 10. Growth rates exhibit a quick convergence to the y ~ 
0.7 with increasing resolution for all Godunov solvers, with the 
TVDLF solver performing notably worse. The dashed lines are 
marking the uncertainty of growth rate measurements. A decline 
is apparent at the highest resolution. 
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Solver 


32x16x16 


64x32x32 


128x64x64 


256x128x128 


512x256x256 


ZEUSuN 


0.58* 


t0.07 


0.65±0.05 


0.70+0.07 


0.68+0.05 


0.62+0.05 


HLLDlin 


0.53* 


tO.07 


0.65±0.07 


0.70+0.05 


0.70+0.03 


0.68+0.05 


HLLDppfy] 


0.60* 


tO.05 


0.70 ±0.06 


0.71+0.05 


0.71+0.04 


0.67+0.04 


TVDLFuN 


0.33* 


t0.05 


0.53*±0.09 


0.65+0.05 


0.70+0.04 


0.67+0.04 


TVDLFppM 


0.50* 


t0.05 


0.63*±0.08 


0.72+0.06 


0.70+0.03 


0.67+0.04 



Table 2. Growth rates for the resolution tests. The asterisks mark the runs which was not able to resolve n=4. 



6. Conclusion 

The inter-code comparison and the convergence test stress the 
importance of treating the Alfven characteristics for MHD 
Riemann solver in MRI simulations. Observing the conse- 
quences implied by numerical dissipation, we obtain a ro- 
bust and accurate Godunov scheme for 3D MHD simulation 
in curvilinear coordina te systems. The a rithme tic average of 
the 4 Godunov fluxes ( iBalsara & Spicer'. '199 9|) is not consis- 
tent w ith the upwind Godunov scheme (Londrillo & Del Zannal 
I2OOOI) an d leads to numerical insta bility. Beside the loop advec- 
tion test (iGardiner & Stonel l2008h . the linear MRI evolution in 
global disks is a crucial test for a consistent MHD method. Our 
3D MHD simulation with ACT show the unphysical behavior 
of a 'chessboard' instability, mentioned for the HLLD solver 
in iMivoshi & Kusanol (|2008|) . The high-order MHD Riemann 
solvers HLLD and ROE, which contain the Alfven character- 
istic, reproduce the 'chessboard' pattern in the classical ACT 
scheme independent of resolution. 

The upwind EMF reconstruction in combination with the 
accurate and robust HLLD MHD Riemann solver represents a 
good tool for future simulations of 3D MHD accretion disks. 
We reproduce the linear MRI analysis in the inter-code com- 
parison. The Lax-Friedrich and HLL solvers demonstrate much 
higher numerical dissipation in comparison to the finite dif- 
ference scheme ZEUS. The high-order MHD Riemann solver 
HLLD and ROE can compensate this dissipation by using more 
MHD characteristics. The Alfven characteristic appears to play 
an important role for MRI. For the linear MRI evolution both 
HLLD and ROE show the same numerical dissipation as the 
finite difference scheme ZEUS. In addition, the HLLD solver 
shows nearly the exact evolution like ROE, even without re- 
solving the slow magneto-sonic characteristic. In combination 
with the piece-wise parabolic method the HLLD solver has the 
earliest convergence compared to the Lax-Friedrich solver and 
ZEUS. 

7. Outlook 

The Prandtl numbers in global simulations for the high-order 
MHD Riemann solver and physical dissipation are very im- 
portant for MRI turbulence in proto-planetary disks. Proper 
MHD tests shall be done, exploring the numerical dissipation 
of the PLUTO code. iFromang & Papaloizoul (|2007|) measured 
the numerical Prandtl number for ZEUS, w h ich is in the region 
from 2 to 8. Recent work bv ISimon et al.l (l2009l) presented a 
Prandtl number about 2 for the Godunov type ATHENA code. 
Nevertheless, in the proto-planetary disks one can expect rather 
smaller Prandtl numbers, Pm < 1 . The consequences of numeri- 
cal Pm for MRI shall be explored. 
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